llikProfile <- function(model, variable, range, res = 100, parallel = FALSE, ncores = getOption("cl.cores", 1), ...) {
  base <- model

  if (missing(range)) {
    bx <- coef(pbcfit)[deparse(variable)]
    sx <- NULL
    try(sx <- sqrt(vcov(pbcfit)[deparse(variable), deparse(variable)]), silent = TRUE)
    if (!is.null(sx)) {
      range <- c(bx - 3 * sx, bx + 3 * sx)
    } else stop("Range must be specified or there must be coef() and vcov() methods for the model class")
  }

  xval <- seq(from = range[1], to = range[2], length.out = res)

  fun <- function(x, base, variable, ...) {
    mnew <- update(base, substitute(. ~ . - var + offset(k * var), list(var = variable, k = x)))
    logLik(mnew, ...)
  }
  env <- environment()

  if (parallel) {
    require(parallel)
    if (.Platform$OS.type == "windows") {
      n <- names(sessionInfo()$otherPkgs)
      cat(n, fill = TRUE)
      cl <- makeCluster(ncores)
      clusterExport(cl, c("base", "xval", "variable", "fun", "n"), envir = env)
      clusterEvalQ(cl, {
        sapply(n, require, character.only = TRUE)
      })
      llik <- parSapply(cl = cl, X = xval, FUN = fun, base = base, variable = variable)
      stopCluster(cl)
    } else {
      llik <- mclapply(X = xval, FUN = fun, base = base, variable = variable, ...)
    }
  } else {
    llik <- sapply(X = xval, FUN = fun, base = base, variable = variable, ...)
  }
  out <- list(llik = unlist(llik), xval = xval, optimalLL = logLik(base))

  plot(x = out$xval, y = out$llik, type = "l", main = paste("Profile of Log Likelihood for", variable),
    xlab = "Constant Coefficient value", ylab = "Log Likelihood")
  abline(h = out$optimalLL - (3.84/2), col = "blue", lwd = 2)
  return(invisible(out))
}

########### S3 methods for extracting log likelihood
## S3 method for objects of class coxme from package 'coxme'
logLik.coxme <- function(object, type = c("Integrated", "Penalized")) {
  type <- match.arg(type)
  LL <- object$loglik[type]
  return(LL)
}

## S3 method for objects of class coxph from package 'survival'
logLik.coxph <- function(object, type = c("Model", "Null")) {
  type <- match.arg(type)
  type <- switch(type, "Model" = 2, "Null" = 1)
  LL <- object$loglik[type]
  return(LL)
}

if (FALSE) {
require(survival)
pbcfit <- coxph(Surv(time, status==2) ~ age + edema + log(bili) + log(albumin) + log(protime), data = pbc)

system.time(ll1 <- llikProfile(model = pbcfit, variable = quote(log(protime)),
  range = c(0, 4), res = 50, FALSE))
system.time(ll4 <- llikProfile(model = pbcfit, variable = quote(log(protime)),
  range = c(0, 4), res = 4000, parallel = TRUE, ncores = 2))

}

